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We present in this article the use of probabilistic background constraints in 
astronomical image deconvolution to approach to a solution as an interval estimate. 
We elaborate our objective - the interval estimate of the unknown object from 
observed data and our approach - monte-carlo experiment and analysis of marginal 
distributions of image values. One-dimensional observation and deconvolution us- 
ing proposed approach are simulated. Confidence intervals reveal the uncertainties 
due to the background constraint are calculated and significance levels for sources 
retrieved from restored images are provided. © 2012 Optical Society of America 

OCIS codes: 100.0100, 100.1830. 
1. Introduction 

LA. Image restoration and deconvolution 

An imaging system receives the object intensity distribution O as input while yields observed 
intensity distribution I as its output. In most cases (more specifically, if the imaging system is ap- 
proximately linear and shift-invariant, and the output data is contaminated by independent additive 
noise) such a system can be described exactly by the following equation: 

Xoo 
0(xi)P(x - xOdxi + N(x) 
(1) 

= (0*P) (x) + N(x) 
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where x is a spatial coordinate variable, P is the point spread function (PSF), N is the additive 
noise and the asterisk * denotes the convolution operation HI. 

The problem of reversing the effects of convolution on the observed data /, namely, restore the 
object intensity O, with an observation /, as well as a determinate PSF P is called deconvolution 
problem, (otherwise if the PSF is not determined we call the problem blind deconvolution problem 
[fJl, which is not discussed in this article). 

l.B. Deconvolution methods and constraints 

Under suitable conditions the Fourier transform of a convolution is the point-wise product of 
Fourier transforms, stated by the convolution theorem. According to this direct restorations through 
Fourier transforms are widely used because the thoughts are straightforward and the restorations 
are non-iterative which usually guarantees the efficiency J2]-[5]|. The presence of noise in the convo- 
lution equation [T| gives the problem a statistics aspect. Therefore an estimate of the object intensity 
is also a solution of the deconvolution problem. If the object O can be approximated by a simpler 
function with a few parameters, parametric fits are the best methods [6j. For example the CLEAN 
algorithm, which is designed for radio imaging [7], approximates the object by a series of point 
sources, each with a pair of parameters i.e. the location and flux, and estimates the parameters 
through iterations. Otherwise the deconvolution problem must be modeled non-parametrically. In 
this case different statistics are used to approach an estimate of the object intensity, or a restoration. 
Landweber method minimizes the difference between the real observation / and the convolution 
of the PSF P and an estimate of the object O uniformly ||8l. It implies a Gaussian model since 
this method gives a least-squares estimate. With a Poisson model where the observation / fol- 
lows a Poisson distribution and the major contribution to the noise N comes from photon noise, 
the Richardson-Lucy algorithm computes the maximum likelihood (ML) estimates of the object 
intensity O iteratively flUEEOl. While in Bayesian statistics the maximum a posteriori (MAP) al- 
gorithm maximizes the conditional probability density of the object intensity to achieve a MAP 
estimate [11]. However, methods introduced above do not produce uncertainty, but only a single 
restoration |[T2l . The expectation through Markov Chain Monte Carlo (EMC2) method as an ex- 
tension of the Richardson-Lucy algorithm uses the same likelihood of the statistical model but 
sampling -based fitting method to produce uncertainty information 021. 

Deconvolution problem is known as an ill-posed problem [13]. Constraints are vital for a desired 
solution of the problem - an image represents the object with acceptable errors. Tighter constraints 
prevent underfitting the observation while looser constraints prevent overfitting it BH. The use of 
constraints is termed regularization ||TTllT4l [T3Tl . Commonly used constraints include positivity, 
backgrounds, lower bounds (i.e. backgrounds) and upper bounds |[T6l . band-limited constraints in 
Fourier domain (as various window functions or filters e.g. Wiener filter), imaging entropy, etc. 
Positivity constraint is implied in Richardson-Lucy algorithm as well as in the MAP algorithm. It 
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helps Richardson-Lucy algorithm avoid amplifying noise, which is an important drawback of this 
algorithm ifTTI . In this article we will focus on background constraints in MAP deconvolution. 

1. C. About this article 

The concept of background constraint is straightforward - the background of an image implies the 
lower bounds of intensities of pixels on the image. In general cases the background is estimated 
before it is used as a constraint. 

On one hand the performance of deconvolution is substantially affected by the estimation of the 
background lfT8ll . It is worth imposing prior constraints very cautiously because a regularization - 
the use of constraints - makes a compromise between ill-posedness and the accuracy of the decon- 
volution problem |[T9l . More restriction brings not only greater stability but also greater chance to 
eliminate correct solutions to the image reconstruction ||6l. For instance, enforcing a background 
constraint higher than the real one could cost a faint source significant loss or even lose the faint 
source, while providing a lower background constraint could make unpredictable artificial struc- 
tures arise from noise, as demonstrated in the simulation section of this article. 

On the other hand deconvolution methods introduced above are degraded because of the lack 
of uncertainty information in their results |fT2|. In other words, they produce point estimate of the 
unknown object such as a least-squares estimate, an ML estimate, or an MAP estimate instead of an 
interval estimate, which is however often interesting in astromony and astrophysics. The presence 
of an interval estimate, i.e., an image as the "best estimate" of the object with its uncertainty, will 
help discover faint sources, measure significance of restored structures, and retrieve rich statistics 
about interesting sources e.g. the fluxes and positions as well as the corresponding confidence 
intervals. 

We suggest using probabilistic background constraints in deconvolution to take the uncertainty 
of the constraint into account rather than to indulge the uncertainty of it. Moreover, it is a crude 
yet refinable approach to the interval estimate of the deconvolution problem. 

We elaborate the use of probabilistic background constraints in Sect. 2. Next a simulated exam- 
ple is presented in Sect. |3]and is discussed in Sect. 4, where we demonstrate how the approach 
we proposed help discover faint sources, measure significance of sources and fetch interesting 
statistics. Finally the conclusions are given in Sect. 5. 

2. Interval estimate with probabilistic background constraints 

2.A. Uncertainties and interval estimate in deconvolution 

Many factors bring uncertainties to the deconvolution result in astronomical image restoration: 
uncertainty of the PSF, fluctuation of photon arrival i.e. the photon noise, the detector noise (e.g., 
as in the case of typical CCD detectors, the readout noise, the dark current, cosmic rays, uniformity 
of response, and bad pixels, etc.) Il20ll and so on. So it makes sense to model the major noise 
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and employ statistical methods in deconvolution. However the selections between different noise 
models and different methods add more uncertainties. Finally even with a proper noise model and a 
suitable deconvolution method there usually are still several adjustable parameters such as stopping 
conditions in iterative methods, or various prior constraints. Therefore uncertainties are inevitable 
in deconvolution. We target on an interval estimate of the unknown object in astronomical image 
deconvolution. Thus uncertainties are limited in confidence intervals. 

2.B. Approach to an interval estimate: monte-carlo experiment and marginal distributions 

In general cases monte-carlo experiments are feasible to calculate an approximate distribution of a 
given random variable. What we want to calculate is the distribution of possible images given the 
observed data and uncertain background constraints. And from such distribution we can find some 
"best estimate" of the image as well as the corresponding confidence interval of it. When taking 
an image as a multivariate random variable, which groups image value on each pixel together as 
individual variables, the image space is so huge, hence most of the draws sampled from this space 
are not likely to permit the given observed data, thus contributions to the distribution from them are 
quite so small, so monte-carlo experiment sampling directly from the space is largely impractical. 

Instead of directly from the image space, we generate draws from solutions of statistical meth- 
ods, e.g., MAP estimates. For each draw we calculate its probability and use this probability to 
weigh its occurrence. On one hand, in view of the reduction of the image space, the extra calcu- 
lations of probabilities have the merits. On the other hand, given the observed data, the solution 
converges towards the MAP estimate or "best estimate" of other statistics among all possible es- 
timates permitted by the enforced constraints. More specifically, the solution usually converges 
towards the mode, given a specific probability density function of each estimate satisfying the con- 
straints, i.e., estimates with similar probability densities are close to each other. However estimates 
with similar but different background constraints are not that close to each other. Thus we believe 
the distribution of these solutions summarizes the distribution of the population (i.e., the set of all 
possible estimates permitted by given constraints) well enough. 

The "goodness" of an estimate calculated with statistical methods such as the likelihood function 
or the posterior probability is related to joint probability density function of the estimate, i.e., 
the image, but not directly related to image values on individual pixels. However we usually pay 
attention to distributions of image values of interesting structures or even specific pixels rather than 
the image as a whole. So the interval estimate we need is not for the whole image but for image 
values of individual pixels. To achieve that we calculate marginal distributions of image values on 
the pixels. 



4 



2.C. Probabilistic background constraints in de convolution 

In astronomical image deconvolution the background arises from both the remote astrophysical 
components e.g. far away galaxies, diffuse sources, etc., and the local instrumental factors. Com- 
pared with PSF, detector noise, and other uncertainties sources the background is more difficult to 
be determined precisely through prior calibrations in a single observation. This article focuses on 
the background constraints only. Here we will not discuss uncertainties from PSF, noise, or other 
aspects. 

While we estimate the background from the observed data we can also estimate the confidence 
interval of the estimated background. The estimation can be achieved by following steps. First, ex- 
tract the observed data dominated by background luminosities. We call it background data. Second, 
according to the background data itself or prior knowledge select a proper background model e.g. 
a spatial function with several adjustable parameters. Finally, calculate the "best-fit" parameters as 
well as their confidence intervals with the background data. 

For example, If the object consists of sparse sources only and the background has only large- 
scale spatial structures, the background can be approximated with a simple function or even a con- 
stant soon after the sources have been removed from the observed data. Clipping methods H18U21I 
and CLEAN algorithm |[T6ll are both workable to extract the background data. We pick CLEAN 
algorithm and assume the image values on the residual map Ir(x), which serves as the extracted 
background data, follows a normal distribution. As a result we can use least-squares method to 
estimate the background parameters as well as their confidence intervals. To make the example as 
simple as possible we assume further that the background does not have spatial structures at all 
so we can use a constant B to model the background. Thus the average of the background data 
Ir(x) is exactly the "best-fit" background, i.e. B = Ir(x) (since the convolution of a constant and a 
normalized function is the constant itself, i.e., (B * P)(x) = B) and the variance of the background 



data Var[I R (x)] = [Ir(x) - Ir(x)) reflects its uncertainty. Given the central limit theorem the es- 
timated background follows a normal distribution approximately, and the mean of the estimated 
background is also B while its standard deviation is related to the variance of the background data 
as erg = (Var[I R (x)]/N) l/2 given the background data consists of N observed count rates. 

In iterative methods such as Richardson-Lucy algorithm the background constraint is enforced 

as: 



image corrected by the background constraint lfTTl[T8l . B usually is the estimated background. 
We suggest using a probabilistic background constraint B k sampled from the distribution of the 
background estimate through a series of iterations, namely, performing the k-th deconvolution with 





otherwise 




(2) 



where 0(x) is an estimate of the unknown object thus an image at a iteration step, 0'{x) is the 
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Fig. 1. Simulated object 0(x). Sources: 1000 photons/ s at 200 pixel, 
4000 photons Is at 400 pixel, 10000 photons /s at 600 pixel, and 5 000 photons /s 
at 800 pixel. The intensity of the background is 1 000 photons js. 



the background constraint B k . This is the monte-carlo experiment we suggest in Sect. |2.B| 

3. Simulation and results 

3. A. Configuration of simulation 
3.A.I. Object, PSF and observed data 

The simulated object consists of several point sources and the background are shown in Fig. [TJ 
The PSF is simulated with a normalized Gaussian bell-curve with FWHM (full width at half 

maximum) of 46 pixels (as Fig. [2]). 

The observed data (e.g. Fig. [3]) is simulated by: 

1. convolving 0(x) with P(x) (let I(x) = (O * P)(x)), 
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Fig. 2. Simulated PSF P(x). 

2. adding photon noise to I(x) (we use normal fluctuation N(0, I(x)) to approximate the photon 
noise with Poisson distributions, Pois( y/(x)).), 

3. adding additive normal noise N a b s (N b S ~ N(0, 25)) to I(x) to simulate other noise sources, 
and 

4. rounding each observed data value to its nearest integer to simulate the round-off error. 

The simulated object as well as the PSF and the observed data are discretized into 1 024 x 1 pixels 
vectors. Thus the coordinate variable x of either O(x), P(x) or I(x) is discretized into a dimension- 
less integer variable in the sense that it denotes an arbitrary pixel. 
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Observed data 




Fig. 3. Simulated observed data. 



3.A.2. Background estimation 

CLEAN algorithm is employed as introduced in Sect. l2.C[ The residual map, i.e. the background 
data is extracted as shown in Fig. @J The arithmetic mean of the data Ir(x) = 1 002.1 counts/ s 

r il/2 

while its standard deviation o- Ir(x) « (I R (x) - Ir(x)) 2 =31.5 counts/ s. 
From the background data I R (x) (Fig. ID) we obtain: 

• the "best estimate" of the background intensity 5= 1 002.1 counts /s, and 

• the estimate of the background is approximately distributed as N(\ 002.1, 1.2). 
3.A.3. MAP deconvolution 

We employ MAP method to solve the deconvolution problem. In this configuration we find within 
2 500 iterations the maximization process converges well enough to provide a stable estimate, i.e. 
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Extracted background data 




Fig. 4. Background data Ir(x). 



Omap, as shown in Fig. [5] 

We also investigate the convergence of this method by observe the mean square error as well as 
the posterior probability density going with the iterations. The mean square error in the k-th step 
is defined as 

^ = pW - O ik M ^{x)f 

and it measures the difference between two adjacent estimates, say, the k-th estimate and its previ- 
ous one. The posterior probability density is defined as |[TT| 

p(o y MAP \i) = n, — (3) 

if assume a uniform prior probability of the unknown object. It is convenient calculating the log- 
arithm to base 10 of the probability density in numerically since even that of the most probable 
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Fig. 5. Object estimate 6 map- Deconvolving with 2 500 MAP iterations, using B = 
I R (x) as the background constraint. 

estimate is quite small. See Fig. [6] and Fig. [7] for the mean square error and posterior probability 
density respectively. 

The posterior probability density shown in Equation [3] is evaluated approximately. 

3.B. Monte-carlo experiment and results 
3.B.I. Monte-carlo experiment 

As suggest in Sect. I2.BI we setup a monte-carlo experiment to calculate marginal distributions of 
image values on individual pixels. The experiment consists of following steps: 

1 . Generate a random draw B k from N{\ 002. 1 , 1 .2), which is the distribution of the background 
estimates. 

2. Use B k as the background constraint in deconvolving the observed data I(x). 
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Fig. 6. Convergence of MSE (mean square error). It shows the MSE between two 
adjacent iteration steps. 



3. As a result we get a MAP solution O k {x) and its probability p(O k \l). 

4. Repeat from 1 to 3 for K times, store all K Ok{x) and p{Ou\I). 
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Fig. 7. Convergence of posterior PDF (probability density function). It shows the 
logarithm to base 10 of the posterior probability density function of the estimate in 
each iteration step. 



See the diagram of monte-carlo experiment below for the setup. 

Distribution of background estimates 




/ MAP method O k , p(O k \I) 
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In the diagram / is the observed data, Ou and p{Ok\I) are the image and its probability given the 
observed data in the k-th deconvolution with MAP method. 

3.B.2. Marginal distribution and interval estimate 

First we calculate the histogram for image value on each pixel. Then we calculate the total posterior 
probability density for each interval of each histogram. Finally we scale the bar on each interval 
by the total posterior probability density to update the histograms. See Fig.[8]for examples. Such 
histograms portray the marginal distributions of image values individually. 

Because the background is not restored by deconvolution process but enforced by the constraint, 
it is truncated from each image value in Fig. [8j 

We find the median as well as the central confidence interval (CI) of the marginal distribution of 
the image value for each pixel from its histogram. For example there are 0.05% image values are 
greater than the upper limit of the 99.9% confidence interval and 0.05% image values are less than 
the lower limit of it. See the medians as well as the confidence intervals in Fig.[9l 

3.B.3. Source retrieval and significance level 

We use a program to search for continuous structures in each image restored from the monte- 
carlo experiment while marginal distributions of image values are being calculated. Each structure 
is considered as a source. We take the sum of all the image values of pixels within the range of a 
source as its flux, while we take the barycentric position, weighted by image values, as its position. 
We account sources overlap with one another between different images the same source. So we 
obtain a distribution for the flux of each source and that for its position. From the distribution, 
we can not only calculate the "best estimate" of each flux or position but also the corresponding 
confidence interval (as shown in Fig. [TTI) . 

We perform a monte-carlo experiment of deconvolutions on random background data sets to 
calculate the null distribution of flux in restored images. The random background data sets have 
the same means and variances as the background data extracted from the observed data in Sect. 
I3.A.2I We retrieve fluxes of "sources" arising from "images" restored from the background data. 
Of course these restored "sources" arise from magnified noise since there is no source in the back- 
ground data at all. Finally we calculate the null distribution of fluxes of these artificial "sources" 
we restored from random background data sets. See the distribution in Fig.fTOl 

Vertical lines in Fig.fTOl indicate significance levels of sources restored from observed data. For 
example there are only 5% artificial "sources" arising from the noise have fluxes greater than the 
flux where the cumulative probability is 95%. So the dashed line denoting this probability also 
indicates the 5% (about 2<x) significance level, i.e., 3 counts/ s. The 1% (about 2.6cr), 0.5% (about 
2.8cr), and 0.1% (about 3.3<r) significance levels are 80 counts/ s, 98 counts/ s and 142 counts / s 
respectively. 
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4. Discussion 

4.A. About results of numerical simulation 

In Sect. [3] we calculated the "best estimate" of the unknown object as well as its confidence inter- 
vals as shown in Fig. [9] The performance of deconvolution is affected by the background constraint 
we enforced. Compared with the object (as in Fig.Q]), there are artificial structures in the restored 
images, and the intensities of these structures vary in images restored with different background 
constraints. In a real case we certainly have no idea whether such a structure in an image restored 
from observed data is artificial and due to the underestimate of the background, or a real yet faint 
source. Given the confidence interval, however, we know the lower bounds. So if the lower bound 
of a structure is still positive we know this structure is not due to the underestimate of the back- 
ground at least. For example, as shown in Fig. [51 there are artifacts in the image restored by a single 
run of deconvolution. They are mainly due to the underestimate of the background since their pres- 
ences are sensitive to the background constraints obviously. For bright sources, their confidence 
intervals reveal uncertainties due to the estimate of the background. 

As shown in Fig.[8]the marginal distributions are highly skewed and some of them are bimodal. 
Compared with other measures of central tendency such as mean and mode, we use the median 
as the best estimate. The mode of such distribution converges only with plenty of scores been 
accumulated. The mean tends to be affected by extreme values. Both the mode and the mean can 
be out of the central confidence interval. We believe the median is good enough to measure the 
central tendency and it is insensitive to the sample size. 

As shown in Fig.QTIfhe fluxes of sources are more dispersed than the positions. This is due to 
the limitation as well as the aim of the approach we take to calculate the confidence intervals. Our 
approach is aimed at and limited to the background constraints. So the confidence intervals are 
all about the background constraints. Therefore it is clear that the background constraint mainly 
affects the intensities of restored image values in MAP deconvolution. Also we see in Fig. ITTI 
even 99.9% confidence intervals do not cover all correct fluxes and positions. It is still because 
the intervals only measure the uncertainties brought by the background constraint. The extra errors 
are mainly caused by simulated photon noise and detector noise. Since we perform deconvolutions 
only on a single set of observed data, errors caused by these noises can not be eliminated. 

4.B. Feasibility of our approach in general cases 

In Sect. [2] and and Sect. [3] we focus our attention on background constraint in deconvolution. 
In Sect. [3] we perform numerical simulations employing MAP estimate. However the approach 
suggest in this article is not confined in the case we elaborated or simulated. 

With necessary prior analysis on the distribution of such as the error of PSF, the detector noise, 
etc., we can always perform monte-carlo experiment of deconvolutions with random draw of cor- 
responding input. Then with the results we can calculate marginal distributions of image values on 



14 



all pixels immediately. Although these already reflect the uncertainty brought by the corresponding 
input, since in our approach the results of monte-carlo experiment are stored individually instead 
of merged prematurely, interesting information can still be retrieved. For example we retrieve the 
fluxes and positions from the restored images in Sect. I3.B.31 It is remarkable that the retrieval pro- 
cedure works as a blind parametric estimate - it is not necessary to know where a source is or how 
many sources there are since it is easier to locate sources in restored images. 

We also suggest running deconvolutions on background data to provide for significance levels 
of sources in images restored in observed data, as shown in Fig. [TOl This is about another source 
of artifacts in image restoration - magnified noise in observed data. We take only the noise on 
background in to account, however, artifacts could also arise nearby a bright source due to the 
noise on its flux spreaded by PSF To analyze significance levels of source restored nearby a bright 
source, the bright source should also be put into the data. 

5. Conclusions 

We elaborated and simulated using probabilistic background constraints to build monte-carlo ex- 
periment of deconvolutions. From the results we calculated interval estimates of restored images 
that reveal the uncertainties brought by the background constraints to the deconvolution. The back- 
ground constraints affect mainly the intensities of restored image values. This approach helps dis- 
criminate structures arising from restored images caused by the underestimate of the background 
from others and discover those disappear due to the overestimate of the background. We can also 
use this approach to analyze uncertainties from other aspects thus eliminate errors from other as- 
pects. 

Acknowledgements 

This work was supported by the National Natural Science Foundation of China (NSFC) under 
Grant No. 1 1 173038, and also by Tsinghua University Initiative Scientific Research Program under 
Grant No. 20111081102. 

References 

1. A. Meister, Deconvolution Problems in Nonparametric Statistics (Springer-Verlag, Berlin, 
2009). 

2. G. R. Ayers and J. C. Dainty, "Iterative blind deconvolution method and its applications," 0l 
13, 547-549 (1988). 

3. S. James L. Harris, "Image Evaluation and Restoration," J. Opt. Soc. Am. 56, 569-574 (1966). 

4. S. Zaroubi, Y. Hoffman, K. B. Fisher, and O. Lahav, "Wiener Reconstruction of the Large- 
Scale Structure," Astrophys. J. 449, 446 (1995). 



15 



5. M. G. Lofdahl, "Multiframe deconvolution with space- variant point-spread functions by use 
of inverse filtering and fast Fourier transform," Appl. Opt. 46, 4686-4693 (2007). 

6. R. Puetter, T. Gosnell, and A. Yahil, "DIGITAL IMAGE RECONSTRUCTION: Deblurring 
and Denoising," Annual Review of Astronomy and Astrophysics 43, 139-194 (2005). 

7. J. A. Hogbom, "Aperture Synthesis with a Non-Regular Distribution of Interferometer Base- 
lines," Astronomy and Astrophysics Supplementl5, 417 — h (1974). 

8. L. Landweber, "An Iteration Formula for Fredholm Integral Equations of the First Kind," 
American Journal of Mathematics 73, pp. 615-624 (1951). 

9. W. H. Richardson, "Bayesian-Based Iterative Method of Image Restoration," J. Opt. Soc. Am. 
62, 55-59 (1972). 

10. L. B. Lucy, "An iterative technique for the rectification of observed distributions," The Astro- 
nomical Journal 79 (1974). 

11. J. L. Starck, E. Pantin, and F. Murtagh, "Deconvolution in Astronomy: A Review," The Publi- 
cations of the Astronomical Society of the Pacific 114, 1051-1069 (2002). 

12. D. N. Esch, A. Connors, M. Karovska, and D. A. van Dyk, "An Image Restoration Technique 
with Error Estimates," Astrophys. J. 610, 1213-1227 (2004). 

13. C. Byrne, "Iterative algorithms for deblurring and deconvolution with constraints," Inverse 
Problems 14, 1455 (1998). 

14. A. Lannes, S. Roques, and M. J. Casanove, "Stabilized Reconstruction in Signal and Image 
Processing - 1. Partial Deconvolution and Spectral Extrapolation with Limited Field," Journal 
of Modern Optics 34 (1987). 

15. P. Hansen, "Regularization Tools," Numerical Algorithms 6, 1-35 (1994). 

16. T.-P. Li and M. Wu, "A Direct Restoration Method for Spectral and Image Analysis," Astro- 
physics and Space Science206, 91-102 (1993). 

17. N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J.-C. Olivo-Marin, and J. Zerubia, 
"Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope 
deconvolution," Microscopy Research and Technique 69, 260-266 (2006). 

18. G. M. P. van Kempen and L. J. van Vliet, "Background estimation in nonlinear image restora- 
tion," J. Opt. Soc. Am. A 17, 425-433 (2000). 

19. G. M. Wing, A Primer on Integral Equations of the First Kind: The Problem of Deconvolution 
and Unfolding (Society for Industrial Mathematics, Philadelphia, 1991). 

20. I. S. McLean, Electronic Imaging in Astronomy: Detectors and Instrumentation (Springer, 
2008), 2nd ed. 

21 . J.-L. Starck and F. Murtagh, Astronomical Image and Data Analysis (Springer- Verlag, Berlin, 
2006). 



16 



8000 



Image value distribution at pixel 310 




10 15 20 

rate counts/s 



Image value distribution at pixel 400 




358 



360 



362 364 366 
rate counts/s 



368 



370 



(b) 

Image value distribution at pixel 600 




1020 



1040 



1060 1080 
rate counts/s 



1100 



1120 



(c) 



Fig. 8. Normalized occurrence frequency of image values (backgrounds are trun- 
cated) on x = 310 pixel (top), x = 400 pixel (middle), and x = 600 pixel (bottom). 
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Image with 99.9% confidence intervals 
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Fig. 9. Medians of marginal distributions of image values for all the pixels as well 
as their 99.9% (about 33a) confidence intervals in images. Medians, upper limits 
and lower limits are indicated by solid line, dashes, and dots respectively. Data 
series are shifted, smoothed and plotted with logarithmic scaled y-axis to highlight 
differences between different series especially on faint structures. 
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Fig. 10. Null distribution of flux restored from random background data sets with 
the same means and variances as the background data extracted from the observed 
data. The dashed line with a circle, the dotted line with a cross, dash-and-dotted line 
with a diamond, and the dashed line with a square denote the 95%, 99%, 99.5%, 
and 99.9% cumulative probabilities respectively. 
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Fig. 1 1 . (Color online) Fluxes (y-coordinates) and positions (^-coordinates) of re- 
trieved sources, 99.9% confidence intervals of both fluxes and positions, as well as 
the 0.01% significance level for fluxes of retrieved sources. Blue crosses indicate 
the fluxes and of sources retrieved from monte-carlo experiment using probabilis- 
tic background constraints, while error bars along jc-axis indicate central 99.9% 
confidence intervals of positions and error bars along v-axis indicate central 99.9% 
confidence intervals of fluxes. Green solid lines with circles indicate sources in the 
simulated object. Red dashed lines with triangles indicate sources retrieved from 
a single run of MAP deconvolution. Black dotted line indicates the 0.01% signifi- 
cance level. Here the significance level is 142 counts I s. 
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